function [rho,rhou,E] = LaxFriedrichs4(n,v)
%用LaxFridrichs格式计算Euler方程的初边值问题4
%参数说明：
%   输入:
%   n：单位长度上的网格数量。
%   v:\tau / h
%   输出：
%   rho,rhou,E:相应变量在区域[-5,5]*[0,T]上的值
    T= 1.8;
    
    h = 10/n;
    tau = h*v;
    
    step = floor(T/tau);
    %初始化
    rho = zeros(n+1,2);
    rhou = zeros(n+1,2);
    E = zeros(n+1,2);
    
    %初始条件
    a = floor(0.1*n);

    rho(1:a,1)=3.857143;
    rhou(1:a,1)=10.14185;
    E(1:a,1)= 39.16666;
    
    rho(a+1:n+1,1)=1+0.2*sin((a+1:n+1)*10/n-5-10/n);
    rhou(a+1:n+1,1)=0;
    E(a+1:n+1,1)= 2.5;
    
    frho=zeros(1,n+1);
    frhou=zeros(1,n+1);
    fE=zeros(1,n+1);
    %迭代
    ti  = 2;
    t = 1;
    for s = 2:step+1
        x = ti;
        ti = t;
        t = x;
        %计算f_j
        for i = 1:n+1
            rhou2 =(rhou(i,ti)^2/rho(i,ti));
            p = 0.4*(E(i,ti)-rhou2/2);
            frho(i)=rhou(i,ti);
            frhou(i)=rhou2+p;
            fE(i)=rhou(i,ti)/rho(i,ti)*(E(i,ti)+p);
        end
        
        for i=2:n
            rho(i,t)=(rho(i-1,ti)+rho(i+1,ti))/2-v/2*(frho(i+1)-frho(i-1));
            rhou(i,t)=(rhou(i-1,ti)+rhou(i+1,ti))/2-v/2*(frhou(i+1)-frhou(i-1));
            E(i,t)=(E(i-1,ti)+E(i+1,ti))/2-v/2*(fE(i+1)-fE(i-1));
        end
        rho(1,t)=rho(2,t);
        rho(n+1,t)=rho(n,t);
        rhou(1,t)=rhou(2,t);
        rhou(n+1,t)=rhou(n,t);
        E(1,t)=E(2,t);
        E(n+1,t)=E(n,t);
    end
    rho=rho(:,t);
    rhou=rhou(:,t);
    E=E(:,t);
    ax = 0:1/n:1;
    plot(ax,rho);
    title(['Question 4 \rho at time 1.8']);
    h = gcf;
    hgexport(h,['report/LF/q4.eps']);
            
end